#### setting environment ####
# install.packages("estimatr")
require(estimatr)
windowsFonts(Helvetica = windowsFont("Helvetica"))

data <- read.csv("replication_data.csv")
additional.data <- read.csv("additional_data.csv")

#### descriptive statistics (Appendix B.2) ####
# number of respondents in each group (Table A1)
table(data$Z, 2 - data$W)
table(2 - data$W)
table(data$Z)

#### balance tests (Appendix C.1) ####
# function to conduct the F-test for equivalence
equivalence.F.test <- function (data,  # data frame containing variables
                                variable,  # variable to be tested
                                group,  # group variable
                                alpha,  # significance level
                                epsilon,  # equivalence limits
                                digit  # number of decimal places
                                ) {
  variable.comp <- na.omit(data[, variable])
  group.comp <- data[! is.na(data[, variable]), group]
  N <- length(variable.comp)
  k <- length(unique(group.comp))
  n <- table(group.comp)
  n.bar <- mean(n)
  grand.mean <- mean(variable.comp)
  group.mean <- tapply(variable.comp, group.comp, mean)
  SS.between <- 0
  SS.within <- rep(NA, k)
  for (i in 1:k) {
    SS.between <- SS.between + 
      (n[i] / n.bar) * (group.mean[i] - grand.mean) ^ 2
    SS.within[i] <- sum((variable.comp[group.comp == i] - group.mean[i]) ^ 2)
  }
  psi.sq.hat <- SS.between / 
    (sum(SS.within) / (N - k))  # Eq.(7.9) in Welleck (2010)
  noncentrality.par <- n.bar * epsilon ^ 2
  threshold <- ((k - 1) / n.bar) * qf(alpha, k - 1, N - k, noncentrality.par)
  test.result <- psi.sq.hat < threshold  # Eq.(7.10) in Welleck (2010)
  result <- data.frame(matrix(round(group.mean, 2), 1, k), test.result)
  rownames(result) <- variable
  result
}

# conduct the F-test for equivalence
equivalence.F.test.result <- 
  rbind(equivalence.F.test(data, "gender", "group", 0.05, 0.25), 
        equivalence.F.test(data, "age", "group", 0.05, 0.25), 
        equivalence.F.test(data, "low.edu", "group", 0.05, 0.25), 
        equivalence.F.test(data, "middle.edu", "group", 0.05, 0.25), 
        equivalence.F.test(data, "high.edu", "group", 0.05, 0.25), 
        equivalence.F.test(data, "income", "group", 0.05, 0.25), 
        equivalence.F.test(data, "knowledge", "group", 0.05, 0.25), 
        equivalence.F.test(data, "government", "group", 0.05, 0.25), 
        equivalence.F.test(data, "opposition", "group", 0.05, 0.25), 
        equivalence.F.test(data, "independent", "group", 0.05, 0.25), 
        equivalence.F.test(data, "ideology", "group", 0.05, 0.25))

# results of the F-test for equivalence (Table A2)
equivalence.F.test.result

#### factual manupulation check (Appendix C.2) ####
# proportion of correct answers in the high-skill condition
round(mean(data$MC.correct[data$W == 1]), 3)
# proportion of correct answers in the low-skill condition
round(mean(data$MC.correct[data$W == 0]), 3)
# proportion of DK in the high-skill condition
round(mean((data$MC == 5)[data$W == 1]), 3)
# proportion of DK in the low-skill condition
round(mean((data$MC == 5)[data$W == 0]), 3)

#### main analysis (main text and Appendix F) ####
## only NMA
lm.skill <- lm_robust(Y ~ W, data = data, 
                      subset = Z.0 == 1)
# Table A3, top left
round(summary(lm.skill)$coefficients, 3)
lm.skill$nobs

# approval rate for a low-skilled immigrant
lm.low.skill <- lm_robust(Y ~ 1, data = data, 
                          subset = Z.0 == 1 & W == 0)
round(summary(lm.low.skill)$coefficients, 3)
# approval rate for a high-skilled immigrant
lm.high.skill <- lm_robust(Y ~ 1, data = data, 
                           subset = Z.0 == 1 & W == 1)
round(summary(lm.high.skill)$coefficients, 3)

## NMA and economic MMA
# linear regression
lm.economic <- lm_robust(Y ~ W * Z.0, data = data, 
                         subset = (Z.0 == 1 | Z.1 == 1))
# Table A3, top right
round(summary(lm.economic)$coefficients, 3)
lm.economic$nobs

## NMA and welfare MMA
lm.welfare <- lm_robust(Y ~ W * Z.0, data = data, 
                        subset = (Z.0 == 1 | Z.2 == 1))
# Table A3, bottom left
round(summary(lm.welfare)$coefficients, 3)
lm.welfare$nobs

## NMA and criminal MMA
lm.criminal <- lm_robust(Y ~ W * Z.0, data = data, 
                         subset = (Z.0 == 1 | Z.3 == 1))
# Table A3, bottom right
round(summary(lm.criminal)$coefficients, 3)
lm.criminal$nobs

## Figure 1
png("Fugire_1.png", 
    width = 6, height = 2.5, units = "in", pointsize = 10, res = 1500)
layout(matrix(1:2, 1, 2))
par(mar = c(3.5, 0, 2.5, 0))
plot(NULL, NULL, type = "n", bty = "n", 
     xlim = c(-0.3, 0.25), ylim = c(0, 4.5), 
     xlab = "", ylab = "", xaxt = "n", yaxt = "n")
mtext("ATE and ACDEs", line = 1, at = 0.05, cex = 1.2, font = 2)
abline(v = 0, col = "gray")
abline(v = c(-0.1, 0.1, 0.2), lty = 3, col = "gray")
points(lm.skill$coefficients[2], 4, pch = 19)
segments(lm.skill$conf.low[2], 4, lm.skill$conf.high[2], 4)
points(lm.economic$coefficients[2], 2.5, pch = 19)
segments(lm.economic$conf.low[2], 2.5, lm.economic$conf.high[2], 2.5)
points(lm.welfare$coefficients[2], 1.5, pch = 19)
segments(lm.welfare$conf.low[2], 1.5, lm.welfare$conf.high[2], 1.5)
points(lm.criminal$coefficients[2], 0.5, pch = 19)
segments(lm.criminal$conf.low[2], 0.5, lm.criminal$conf.high[2], 0.5)
text(-0.3, 4, "ATE", pos = 4)
text(-0.3, 3, "ACDE", pos = 4)
text(-0.25, 2.5, "Economic", pos = 4)
text(-0.25, 1.5, "Welfare", pos = 4)
text(-0.25, 0.5, "Criminal", pos = 4)
axis(1, at = seq(-0.1, 0.2, 0.1))
mtext("Effects", side = 1, at = 0.05, line = 2.5)
plot(NULL, NULL, type = "n", bty = "n", 
     xlim = c(-0.35, 0.2), ylim = c(0, 4.5), 
     xlab = "", ylab = "", xaxt = "n", yaxt = "n")
mtext("EEs", line = 1, at = 0.05, cex = 1.2, font = 2)
abline(v = 0, col = "gray")
abline(v = c(-0.1, 0.1, 0.2), lty = 3, col = "gray")
points(lm.economic$coefficients[4], 2.5, pch = 19)
segments(lm.economic$conf.low[4], 2.5, lm.economic$conf.high[4], 2.5)
points(lm.welfare$coefficients[4], 1.5, pch = 19)
segments(lm.welfare$conf.low[4], 1.5, lm.welfare$conf.high[4], 1.5)
points(lm.criminal$coefficients[4], 0.5, pch = 19)
segments(lm.criminal$conf.low[4], 0.5, lm.criminal$conf.high[4], 0.5)
text(-0.35, 3, "EE", pos = 4)
text(-0.3, 2.5, "Economic", pos = 4)
text(-0.3, 1.5, "Welfare", pos = 4)
text(-0.3, 0.5, "Criminal", pos = 4)
axis(1, at = seq(-0.1, 0.2, 0.1))
mtext("Effects", side = 1, at = 0.05, line = 2.5)
dev.off()

#### additional experiment (Appendix H) ####
# number of respondents
nrow(additional.data)

## average across treatment groups
barplot.data <- array(NA, c(2, 6, 3))
variable.names <- paste("Y", 1:6, sep = ".")
for (i in 1:6) {
  temp.data <- data.frame(Y = additional.data[, variable.names[i]], 
                          W = additional.data$W)
  high.skill <- lm_robust(Y ~ 1, data = temp.data, subset = W == 1)
  barplot.data[1, i, 1] <- high.skill$coefficients
  barplot.data[1, i, 2] <- high.skill$conf.low
  barplot.data[1, i, 3] <- high.skill$conf.high
  low.skill <- lm_robust(Y ~ 1, data = temp.data, subset = W == 0)
  barplot.data[2, i, 1] <- low.skill$coefficients
  barplot.data[2, i, 2] <- low.skill$conf.low
  barplot.data[2, i, 3] <- low.skill$conf.high
}

## Figure A1
png("Fugire_A1.png", 
    width = 6, height = 4, units = "in", pointsize = 8, res = 1500)
par(mar = c(3, 3.5, 0.5, 0.5), lwd = 0.5)
plot(NULL, NULL, type = "n", bty = "n", 
     xlim = c(0.5, 18.5), ylim = c(0, 4), 
     xlab = "", ylab = "", xaxt = "n", yaxt = "n")
abline(h = seq(0, 4, 1), lty = 3, col = "gray")
barplot.result <- barplot(barplot.data[, , 1], beside = TRUE, 
                          col = c("gray50", "white"), ylim = c(0, 4), 
                          xlab = "", ylab = "", xaxt = "n", yaxt = "n", 
                          add = TRUE)
segments(barplot.result[1, ], barplot.data[1, , 2], 
         barplot.result[1, ], barplot.data[1, , 3])
segments(barplot.result[2, ], barplot.data[2, , 2], 
         barplot.result[2, ], barplot.data[2, , 3])
mtext(rep("H", 6), side = 1, at = barplot.result[1, ], line = -0.5)
mtext(rep("L", 6), side = 1, at = barplot.result[2, ], line = -0.5)
mtext(c("Financial\nanalyist", "Business\nconsultant", "Think tank\nresearcher", 
        "Convenience\nstore clerk", "Construction\nworker", "Janitor\n"), 
      side = 1, at = colMeans(barplot.result), line = 2)
axis(2, lwd = 0.5)
mtext("Perceived likelihood", side = 2, line = 2.5)
dev.off()

## statistical significance
lm.financial.analyist <- lm_robust(Y.1 ~ W, data = additional.data)
round(summary(lm.financial.analyist)$coefficients, 3)

lm.business.consultant <- lm_robust(Y.2 ~ W, data = additional.data)
round(summary(lm.business.consultant)$coefficients, 3)

lm.think.tank.researcher <- lm_robust(Y.3 ~ W, data = additional.data)
round(summary(lm.think.tank.researcher)$coefficients, 3)

lm.store.clerk <- lm_robust(Y.4 ~ W, data = additional.data)
round(summary(lm.store.clerk)$coefficients, 3)

lm.construction.worker <- lm_robust(Y.5 ~ W, data = additional.data)
round(summary(lm.construction.worker)$coefficients, 3)

lm.janitor <- lm_robust(Y.6 ~ W, data = additional.data)
round(summary(lm.janitor)$coefficients, 3)